#################################################################################
# Representation, Competing Principals, and Waffling on Bills                   #
# in U.S. Legislatures                                                          #
#                                                                               #
# Justin Kirkland and Jeff Harden                                               #
#                                                                               # 
# Empirical replication file                                                    #
# Last update: 6/25/15                                                          # 
#################################################################################
### Packages ###
library(arm)
library(MASS)
library(foreign)
library(ggplot2)

### Data ###
## State legislative data ##
wr <- read.csv("StateLegislatorWaffles20113.csv")
id <- read.csv("StateAssemblyDistrictIDmaster.csv")

# Fill in minority party median for states without minority party leaders (http://www.ncsl.org/legislators-staff/legislators/legislative-leaders/2013-state-legislative-leaders.aspx) #
# id$RepLead_ID[which(id$State == "MS")] <- 0.956
# id$DemLead_ID[which(id$State == "LA")] <- -0.03

# Merge data #
wr.id <- merge(wr, id, by = c("state", "dists"))

## New variables ##
# Party ID #
wr.id$Democrat <- ifelse(wr.id$party == "Democratic" | wr.id$party == "Democratic-Farmer-Labor" | wr.id$party == "Democrat" | wr.id$party == "Progressive", 1, NA)
wr.id$Democrat <- ifelse(wr.id$party == "Republican" | wr.id$party == "Republican/Democratic", 0, wr.id$Democrat)

# Majority party membership ID #
wr.id$inparty <- ifelse(wr.id$majority == "Democrat" & wr.id$Democrat == 1, 1, 0)
wr.id$inparty <- ifelse(wr.id$majority == "Republican" & wr.id$Democrat == 0, 1, wr.id$inparty)

# Extremity measure #
wr.id$extremity <- ifelse(wr.id$Democrat == 1, wr.id$mrp_estimate*-1, wr.id$mrp_estimate)
wr.id$extremity.lead <- ifelse(wr.id$Democrat == 1, wr.id$DemLead_ID*-1, wr.id$RepLead_ID)
wr.id$ex.sq <- wr.id$extremity^2

# Germaneness requirements (http://www.ncsl.org/research/about-state-legislatures/germaneness-requirements.aspx) #
wr.id$germane <- ifelse(wr.id$state == "ms" | wr.id$state == "me" | wr.id$state == "nc" | wr.id$state == "vt" | wr.id$state == "ar" | wr.id$state == "ma" | wr.id$state == "oh" | wr.id$state == "ct" | wr.id$state == "nh" | wr.id$state == "ri", 0, 1)

## Congress data ##
waffle112 <- read.dta("waffle112.dta")
waffle112$leg.id <- paste(waffle112$Lastme, waffle112$State, waffle112$District, sep = "")

d <- subset(waffle112, Party == "Democrat")
r <- subset(waffle112, Party == "Republican")

### Functions ###
## Logit probability function w/ simulation ##
logit.pp <- function(b, vcv, nd, m = 1000, re1 = 0, re2 = 0){

sim.b <- mvrnorm(m, b, vcv)

pred <- matrix(NA, nrow = m, ncol = nrow(nd))
for(j in 1:nrow(nd)){ 
for(i in 1:m){
pred[i, j] <- as.matrix(nd)[j, ] %*% as.matrix(sim.b)[i, ] + re1 + re2
}
}

pe <- apply(plogis(pred), 2, mean)
lo <- apply(plogis(pred), 2, quantile, prob = .025)
hi <- apply(plogis(pred), 2, quantile, prob = .975)

return(list(pe = pe, lo = lo, hi = hi))
}

### Table 1 ###
m.maj <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 1), ], family = binomial (link = logit))

m.min <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 0), ], family = binomial (link = logit))

### Table 2 ###
m.maj.ss <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 1 & wr.id$germane == 1), ], family = binomial (link = logit))

m.min.ss <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 0 & wr.id$germane == 1), ], family = binomial (link = logit))

m.maj.nss <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 1 & wr.id$germane == 0), ], family = binomial (link = logit))

m.min.nss <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 0 & wr.id$germane == 0), ], family = binomial (link = logit))

### Congress Table (Appendix) ###
## Predicted logit coefficients from the state legislative models ##
# From e-mail from Boris Shor, 6/18/15:
# "I have Pelosi at -0.737, Boehner at 0.646, and Cantor at 0.727."
# Prediction for Democrats
d.lead <- -.737*-1 # Multiply by -1 to compute the extremity measure
d.b <- fixef(m.min)[2] + fixef(m.min)[4]*d.lead # Coefficient
d.se <- sqrt(vcov(m.min)[2, 2] + d.lead^2*vcov(m.min)[4, 4] + 2*d.lead*vcov(m.min)[2, 4]) # SE
d.lo <- d.b - 1.96*d.se
d.hi <- d.b + 1.96*d.se
# Coefficient = 0.277864; SE = 0.1386458; 95% CI = [0.006118268, 0.5496096]

# Prediction for Republicans
r.lead <- mean(c(.646, .727))
r.b <- fixef(m.maj)[2] + fixef(m.maj)[4]*r.lead # Coefficient
r.se <- sqrt(vcov(m.maj)[2, 2] + r.lead^2*vcov(m.maj)[4, 4] + 2*r.lead*vcov(m.maj)[2, 4]) # SE
r.lo <- r.b - 1.96*r.se
r.hi <- r.b + 1.96*r.se
# Coefficient = 0.2297216; SE = 0.1481736; 95% CI = [-0.0606986, 0.5201417]

## Models (Table A.5 in the appendix) ##
# Democrats #
m.d <- glmer(matrix(c(SponsoredVotes - BillsNotWaffled, BillsNotWaffled), ncol = 2) ~ ideoex + (1|State) + (1|leg.id), data = d, family = binomial (link = logit))
d.b.hat <- fixef(m.d)[2]
d.se.hat <- sqrt(vcov(m.d)[2, 2])
d.lo.hat <- d.b.hat - 1.96*d.se.hat
d.hi.hat <- d.b.hat + 1.96*d.se.hat
# Coefficient = 0.9992662; SE = 0.4222531; 95% CI = [0.1716501, 1.826882]

# Republicans #
m.r <- glmer(matrix(c(SponsoredVotes - BillsNotWaffled, BillsNotWaffled), ncol = 2) ~ ideoex + (1|State) + (1|leg.id), data = r, family = binomial (link = logit))
r.b.hat <- fixef(m.r)[2]
r.se.hat <- sqrt(vcov(m.r)[2, 2])
r.lo.hat <- r.b.hat - 1.96*r.se.hat
r.hi.hat <- r.b.hat + 1.96*r.se.hat
# Coefficient = 0.1352637; SE = 0.7524027; 95% CI = [-1.339446, 1.609973]
       
### Figures ###
## Figure 3 ##
set.seed(6342)
range.ext.maj <- seq(min(m.maj@frame[ , "extremity"]), max(m.maj@frame[ , "extremity"]), length = 10)
range.ext.min <- seq(min(m.min@frame[ , "extremity"]), max(m.min@frame[ , "extremity"]), length = 10)

pp.ideo1 <- logit.pp(fixef(m.maj), vcov(m.maj), data.frame(cons = 1, ext = range.ext.maj, ext.lead = min(m.maj@frame[ , "extremity.lead"]), int = range.ext.maj*min(m.maj@frame[ , "extremity.lead"])), re1 = median(ranef(m.maj)$state[ , 1]), re2 = median(ranef(m.maj)$leg.s[ , 1])) 
                                   
d.ideo1 <- data.frame(xaxis = range.ext.maj, pe = pp.ideo1$pe, lo = pp.ideo1$lo, hi = pp.ideo1$hi)

pdf("predictwaffle1.pdf")

ggplot(d.ideo1, aes(x = xaxis, y = pe)) + 
geom_line(stat = "identity", lwd = 2, color = "gray50") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = .15) +
scale_y_continuous(limits = c(0, .5), breaks = seq(0, .5, .05)) +
scale_x_continuous(limits = c(min(range.ext.maj), max(range.ext.maj)), breaks = quantile(m.maj@frame[ , "extremity"], c(0, .1, .5, .9, 1)), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(axis.text = element_text(size = 20), axis.title.y = element_text(size = 25, vjust 
= 1.5), axis.title.x = element_text(size = 25, vjust = -.1)) + 
geom_rug(data = m.maj@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

pp.ideo2 <- logit.pp(fixef(m.maj), vcov(m.maj), data.frame(cons = 1, ext = range.ext.maj, ext.lead = max(m.maj@frame[ , "extremity.lead"]), int = range.ext.maj*max(m.maj@frame[ , "extremity.lead"])), re1 = median(ranef(m.maj)$state[ , 1]), re2 = median(ranef(m.maj)$leg.s[ , 1])) 
                                   
d.ideo2 <- data.frame(xaxis = range.ext.maj, pe = pp.ideo2$pe, lo = pp.ideo2$lo, hi = pp.ideo2$hi)

pdf("predictwaffle2.pdf")
 
ggplot(d.ideo2, aes(x = xaxis, y = pe)) + 
geom_line(stat = "identity", lwd = 2, color = "gray50") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = .15) +
scale_y_continuous(limits = c(0, .5), breaks = seq(0, .5, .05)) +
scale_x_continuous(limits = c(min(range.ext.maj), max(range.ext.maj)), breaks = quantile(m.maj@frame[ , "extremity"], c(0, .1, .5, .9, 1)), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(axis.text = element_text(size = 20), axis.title.y = element_text(size = 25, vjust 
= 1.5), axis.title.x = element_text(size = 25, vjust = -.1)) + 
geom_rug(data = m.maj@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

pp.ideo3 <- logit.pp(fixef(m.min), vcov(m.min), data.frame(cons = 1, ext = range.ext.min, ext.lead = min(m.min@frame[ , "extremity.lead"]), int = range.ext.min*min(m.min@frame[ , "extremity.lead"])), re1 = median(ranef(m.min)$state[ , 1]), re2 = median(ranef(m.min)$leg.s[ , 1])) 
                                   
d.ideo3 <- data.frame(xaxis = range.ext.min, pe = pp.ideo3$pe, lo = pp.ideo3$lo, hi = pp.ideo3$hi)

pdf("predictwaffle3.pdf")

ggplot(d.ideo3, aes(x = xaxis, y = pe)) + 
geom_line(stat = "identity", lwd = 2, color = "gray50") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = .15) +
scale_y_continuous(limits = c(0, .5), breaks = seq(0, .5, .05)) +
scale_x_continuous(limits = c(min(range.ext.min), max(range.ext.min)), breaks = quantile(m.min@frame[ , "extremity"], c(0, .1, .5, .9, 1)), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(axis.text = element_text(size = 20), axis.title.y = element_text(size = 25, vjust 
= 1.5), axis.title.x = element_text(size = 25, vjust = -.1)) + 
geom_rug(data = m.min@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

pp.ideo4 <- logit.pp(fixef(m.min), vcov(m.min), data.frame(cons = 1, ext = range.ext.min, ext.lead = max(m.min@frame[ , "extremity.lead"]), int = range.ext.min*max(m.min@frame[ , "extremity.lead"])), re1 = median(ranef(m.min)$state[ , 1]), re2 = median(ranef(m.min)$leg.s[ , 1])) 
                                   
d.ideo4 <- data.frame(xaxis = range.ext.min, pe = pp.ideo4$pe, lo = pp.ideo4$lo, hi = pp.ideo4$hi)

pdf("predictwaffle4.pdf")
 
ggplot(d.ideo4, aes(x = xaxis, y = pe)) + 
geom_line(stat = "identity", lwd = 2, color = "gray50") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = .15) +
scale_y_continuous(limits = c(0, .5), breaks = seq(0, .5, .05)) +
scale_x_continuous(limits = c(min(range.ext.min), max(range.ext.min)), breaks = quantile(m.min@frame[ , "extremity"], c(0, .1, .5, .9, 1)), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(axis.text = element_text(size = 20), axis.title.y = element_text(size = 25, vjust 
= 1.5), axis.title.x = element_text(size = 25, vjust = -.1)) + 
geom_rug(data = m.min@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

## Figure 4 ##
range.ext <- c(min(wr.id$extremity.lead, na.rm = TRUE), max(wr.id$extremity.lead, na.rm = TRUE))

# Single subject provisions #
pp.ideo.ss1 <- logit.pp(fixef(m.maj.ss), vcov(m.maj.ss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 0, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 0, na.rm = TRUE)), re1 = median(ranef(m.maj.ss)$state[ , 1]), re2 = median(ranef(m.maj.ss)$leg.s[ , 1])) 
pp.ideo.ss2 <- logit.pp(fixef(m.maj.ss), vcov(m.maj.ss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 1, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 1, na.rm = TRUE)), re1 = median(ranef(m.maj.ss)$state[ , 1]), re2 = median(ranef(m.maj.ss)$leg.s[ , 1])) 

pp.ideo.ss3 <- logit.pp(fixef(m.min.ss), vcov(m.min.ss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 0, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 0, na.rm = TRUE)), re1 = median(ranef(m.min.ss)$state[ , 1]), re2 = median(ranef(m.min.ss)$leg.s[ , 1])) 
pp.ideo.ss4 <- logit.pp(fixef(m.min.ss), vcov(m.min.ss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 1, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 1, na.rm = TRUE)), re1 = median(ranef(m.min.ss)$state[ , 1]), re2 = median(ranef(m.min.ss)$leg.s[ , 1])) 

# No single subject provisions #
pp.ideo.nss1 <- logit.pp(fixef(m.maj.nss), vcov(m.maj.nss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 0, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 0, na.rm = TRUE)), re1 = median(ranef(m.maj.nss)$state[ , 1]), re2 = median(ranef(m.maj.nss)$leg.s[ , 1])) 
pp.ideo.nss2 <- logit.pp(fixef(m.maj.nss), vcov(m.maj.nss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 1, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 1, na.rm = TRUE)), re1 = median(ranef(m.maj.nss)$state[ , 1]), re2 = median(ranef(m.maj.nss)$leg.s[ , 1])) 

pp.ideo.nss3 <- logit.pp(fixef(m.min.nss), vcov(m.min.nss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 0, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 0, na.rm = TRUE)), re1 = median(ranef(m.min.nss)$state[ , 1]), re2 = median(ranef(m.min.nss)$leg.s[ , 1])) 
pp.ideo.nss4 <- logit.pp(fixef(m.min.nss), vcov(m.min.nss), data.frame(cons = 1, ext = quantile(wr.id$extremity, 1, na.rm = TRUE), ext.lead = range.ext, int = range.ext*quantile(wr.id$extremity, 1, na.rm = TRUE)), re1 = median(ranef(m.min.nss)$state[ , 1]), re2 = median(ranef(m.min.nss)$leg.s[ , 1])) 

# Put data together # 
d.ideo.ss <- data.frame(xaxis = c(.9, 1.1, 1.9, 2.1, 2.9, 3.1, 3.9, 4.1), 
  pe = c(pp.ideo.ss1$pe, pp.ideo.ss2$pe, pp.ideo.ss3$pe, pp.ideo.ss4$pe),
  lo = c(pp.ideo.ss1$lo, pp.ideo.ss2$lo, pp.ideo.ss3$lo, pp.ideo.ss4$lo),
  hi = c(pp.ideo.ss1$hi, pp.ideo.ss2$hi, pp.ideo.ss3$hi, pp.ideo.ss4$hi), 
  party = rep(c("Majority Party", "Minority Party"), each = 4),
  lead = rep(c("Moderate Leaders", "Extreme Leaders"), times = 4),
  dist = rep(c("Moderate District", "Moderate District", "Extreme District", "Extreme District"), times = 2))
d.ideo.ss$dist <- factor(d.ideo.ss$dist, levels = c("Moderate District", "Extreme District"))

d.ideo.nss <- data.frame(xaxis = c(.9, 1.1, 1.9, 2.1, 2.9, 3.1, 3.9, 4.1), 
  pe = c(pp.ideo.nss1$pe, pp.ideo.nss2$pe, pp.ideo.nss3$pe, pp.ideo.nss4$pe),
  lo = c(pp.ideo.nss1$lo, pp.ideo.nss2$lo, pp.ideo.nss3$lo, pp.ideo.nss4$lo),
  hi = c(pp.ideo.nss1$hi, pp.ideo.nss2$hi, pp.ideo.nss3$hi, pp.ideo.nss4$hi), 
  party = rep(c("Majority Party", "Minority Party"), each = 4),
  lead = rep(c("Moderate Leaders", "Extreme Leaders"), times = 4),
  dist = rep(c("Moderate District", "Moderate District", "Extreme District", "Extreme District"), times = 2))
d.ideo.nss$dist <- factor(d.ideo.nss$dist, levels = c("Moderate District", "Extreme District"))

pdf("ss-effects.pdf")
  
ggplot(d.ideo.ss, aes(x = xaxis, y = pe, color = lead)) + 
scale_color_manual(values = c("gray40", "gray60")) +
geom_hline(yintercept = 0, lwd = .5) +
geom_vline(xintercept = 2.5, lwd = .5) +
geom_point(lwd = 4, pch = 19) + 
geom_errorbar(aes(ymin = lo, ymax = hi), size = 1, width = 0, position = position_dodge(width = .25)) +
scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, .1)) +
scale_x_continuous(limits = c(.75, 4.25), breaks = 1:4, labels = c("Moderate District", "Extreme District", "Moderate District", "Extreme District")) +
theme(legend.title = element_blank(), legend.position = c(0, 1), legend.justification = c(0, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 12), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 17, vjust = -.25)) + labs(color = "") + 
ylab("Expected Probability of Waffling") + 
xlab("Majority Party                            Minority Party") 

dev.off()

pdf("nss-effects.pdf")
  
ggplot(d.ideo.nss, aes(x = xaxis, y = pe, color = lead)) + 
scale_color_manual(values = c("gray40", "gray60")) +
geom_hline(yintercept = 0, lwd = .5) +
geom_vline(xintercept = 2.5, lwd = .5) +
geom_point(lwd = 4, pch = 19) + 
geom_errorbar(aes(ymin = lo, ymax = hi), size = 1, width = 0, position = position_dodge(width = .25)) +
scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, .1)) +
scale_x_continuous(limits = c(.75, 4.25), breaks = 1:4, labels = c("Moderate District", "Extreme District", "Moderate District", "Extreme District")) +
theme(legend.title = element_blank(), legend.position = c(0, 1), legend.justification = c(0, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 12), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 17, vjust = -.25)) + labs(color = "") + 
ylab("Expected Probability of Waffling") + 
xlab("Majority Party                            Minority Party") 

dev.off()

## Figure 5 ##
# Democrats #
range.d <- seq(min(d$ideoex, na.rm = TRUE), max(d$ideoex, na.rm = TRUE), length = 10)
pp.d <- logit.pp(fixef(m.d), vcov(m.d), data.frame(cons = 1, ext = range.d), re1 = median(ranef(m.d)$State[ , 1])) 

d.d <- data.frame(xaxis = range.d, pe = pp.d$pe, lo = pp.d$lo, hi = pp.d$hi)

pdf("dems112.pdf")

ggplot(d.d, aes(x = xaxis, y = pe)) + 
geom_line(stat = "identity", lwd = 2, color = "gray50") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = .15) +
scale_y_continuous(limits = c(0, .25), breaks = seq(0, .25, .05)) +
scale_x_continuous(limits = c(min(range.d), max(range.d)), breaks = quantile(m.d@frame[ , "ideoex"], c(0, .1, .5, .9, 1)), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(axis.text = element_text(size = 20), axis.title.y = element_text(size = 25, vjust 
= 1.5), axis.title.x = element_text(size = 25, vjust = -.1)) + 
geom_rug(data = m.maj@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

# Republicans #
range.r <- seq(min(r$ideoex, na.rm = TRUE), max(r$ideoex, na.rm = TRUE), length = 10)
pp.r <- logit.pp(fixef(m.r), vcov(m.r), data.frame(cons = 1, ext = range.r), re1 = median(ranef(m.r)$State[ , 1])) 

d.r <- data.frame(xaxis = range.r, pe = pp.r$pe, lo = pp.r$lo, hi = pp.r$hi)

pdf("reps112.pdf")

ggplot(d.r, aes(x = xaxis, y = pe)) + 
geom_line(stat = "identity", lwd = 2, color = "gray50") +
geom_ribbon(aes(ymin = lo, ymax = hi), alpha = .15) +
scale_y_continuous(limits = c(0, .25), breaks = seq(0, .25, .05)) +
scale_x_continuous(limits = c(min(range.r), max(range.r)), breaks = quantile(m.r@frame[ , "ideoex"], c(0, .1, .5, .9, 1)), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(axis.text = element_text(size = 20), axis.title.y = element_text(size = 25, vjust 
= 1.5), axis.title.x = element_text(size = 25, vjust = -.1)) + 
geom_rug(data = m.maj@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

### Appendix ###
## Table A.1 ##
m.maj2 <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
              extremity + ex.sq + 
              extremity.lead +
              extremity*extremity.lead +
              ex.sq*extremity.lead +
              (1|state) +
              (1|leg.s),
              data = wr.id[which(wr.id$inparty == 1), ], family = binomial (link = logit))

m.min2 <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
              extremity + ex.sq + 
              extremity.lead +
              extremity*extremity.lead +
              ex.sq*extremity.lead +              
              (1|state) +
              (1|leg.s),
              data = wr.id[which(wr.id$inparty == 0), ], family = binomial (link = logit))

# Model fit comparison #
BIC(m.maj, m.maj2); BIC(m.min, m.min2)

## Test for curvilinear relationship ##
ex.range.maj <- seq(min(m.maj2@frame$extremity), max(m.maj2@frame$extremity), length = 10)
ex.range.min <- seq(min(m.min2@frame$extremity), max(m.min2@frame$extremity), length = 10)

d.maj.mod <- data.frame(cons = 1, ex = ex.range.maj, ex.sq = ex.range.maj^2, ex.lead = min(m.maj2@frame$extremity.lead), int1 = ex.range.maj*min(m.maj2@frame$extremity.lead), int2 = (ex.range.maj^2)*min(m.maj2@frame$extremity.lead))
d.maj.ext <- data.frame(cons = 1, ex = ex.range.maj, ex.sq = ex.range.maj^2, ex.lead = max(m.maj2@frame$extremity.lead), int1 = ex.range.maj*max(m.maj2@frame$extremity.lead), int2 = (ex.range.maj^2)*max(m.maj2@frame$extremity.lead))

d.min.mod <- data.frame(cons = 1, ex = ex.range.min, ex.sq = ex.range.min^2, ex.lead = min(m.min2@frame$extremity.lead), int1 = ex.range.min*min(m.min2@frame$extremity.lead), int2 = (ex.range.min^2)*min(m.min2@frame$extremity.lead))
d.min.ext <- data.frame(cons = 1, ex = ex.range.min, ex.sq = ex.range.min^2, ex.lead = max(m.min2@frame$extremity.lead), int1 = ex.range.min*max(m.min2@frame$extremity.lead), int2 = (ex.range.min^2)*max(m.min2@frame$extremity.lead))

pp.maj.mod <- logit.pp(fixef(m.maj2), vcov(m.maj2), d.maj.mod)
pp.maj.ext <- logit.pp(fixef(m.maj2), vcov(m.maj2), d.maj.ext)

pp.min.mod <- logit.pp(fixef(m.min2), vcov(m.min2), d.min.mod)
pp.min.ext <- logit.pp(fixef(m.min2), vcov(m.min2), d.min.ext)

pp.mod <- data.frame(xaxis = c(ex.range.maj, ex.range.min), pe = c(pp.maj.mod$pe, pp.min.mod$pe), lo = c(pp.maj.mod$pe - 1.96*pp.maj.mod$se, pp.min.mod$pe - 1.96*pp.min.mod$pe), hi = c(pp.maj.mod$pe + 1.96*pp.maj.mod$se, pp.min.mod$pe + 1.96*pp.min.mod$pe), label = rep(c("Majority Party", "Minority Party"), each = 10))

# Figure A.2 #
pdf("curvilinear-mod.pdf")

ggplot(pp.mod, aes(x = xaxis, y = pe, color = label)) +
geom_line(lwd = 2) +
scale_color_manual(values = c("gray40", "gray60")) +
scale_y_continuous(limits = c(0, .35), breaks = seq(0, .35, .05)) +
scale_x_continuous(limits = c(min(range.ext.maj), max(range.ext.maj)), breaks = quantile(wr.id$extremity, c(0, .1, .5, .9, 1), na.rm = TRUE), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(legend.title = element_blank(), legend.position = c(0, 1), legend.justification = c(0, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 12), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 17, vjust = -.25)) + labs(color = "") + 
geom_rug(data = m.maj@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

pp.ext <- data.frame(xaxis = c(ex.range.maj, ex.range.min), pe = c(pp.maj.ext$pe, pp.min.ext$pe), lo = c(pp.maj.ext$pe - 1.96*pp.maj.ext$se, pp.min.ext$pe - 1.96*pp.min.ext$pe), hi = c(pp.maj.ext$pe + 1.96*pp.maj.ext$se, pp.min.ext$pe + 1.96*pp.min.ext$pe), label = rep(c("Majority Party", "Minority Party"), each = 10))

pdf("curvilinear-ext.pdf")

ggplot(pp.ext, aes(x = xaxis, y = pe, color = label)) +
geom_line(lwd = 2) +
scale_color_manual(values = c("gray40", "gray60")) +
scale_y_continuous(limits = c(0, .35), breaks = seq(0, .35, .05)) +
scale_x_continuous(limits = c(min(range.ext.maj), max(range.ext.maj)), breaks = quantile(wr.id$extremity, c(0, .1, .5, .9, 1), na.rm = TRUE), labels = c("Min.", "10%", "50%", "90%", "Max.")) +
ylab("Expected Probability of Waffling") + xlab("District Ideological Extremity") +
theme(legend.title = element_blank(), legend.position = c(0, 1), legend.justification = c(0, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text = element_text(size = 12), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 17, vjust = -.25)) + labs(color = "") + 
geom_rug(data = m.maj@frame, aes(x = extremity, y = extremity), position = "jitter", sides = "b", color = "gray50") 

dev.off()

## Figure A.3 ##
by(wr.id$Waffle.Rate, wr.id$state, FUN = mean, na.rm = TRUE)
st <- unique(wr.id$state)
wf.d <- numeric(length(st))
wf.r <- numeric(length(st))

for(i in 1:length(st)){
wf.d[i] <- mean(wr.id$Waffle.Rate[which(wr.id$state == st[i] & wr.id$Democrat == 1)], na.rm = TRUE)  
wf.r[i] <- mean(wr.id$Waffle.Rate[which(wr.id$state == st[i] & wr.id$Democrat == 0)], na.rm = TRUE) 
}

wf.d <- wf.d[-1]
wf.r <- wf.r[-1]
st <- st[-1]

d.wf <- data.frame(wr = c(wf.d, wf.r), party = rep(c("Democrats", "Republicans"), each = length(st)), state = rep(toupper(st), times = 2))

pdf("WaffleByState.pdf", width = 12)

ggplot(d.wf, aes(x = state, y = wr, fill = party)) + 
geom_bar(stat = "identity", lwd = 2, position = position_dodge(width = .95)) +
scale_fill_manual(values = c("gray40", "gray60")) +
scale_y_continuous(limits = c(0, .35), breaks = seq(0, .35, .05)) +
ylab("Mean Waffle Rate") + xlab("State") +
theme(legend.title = element_blank(), legend.position = c(0, 1), legend.justification = c(0, 1), legend.direction = "vertical", legend.text = element_text(size = 15), axis.text.y = element_text(size = 12), axis.text.x = element_text(size = 11, angle = 45, hjust = 1), axis.title.y = element_text(size = 18, vjust = 1.5), axis.title.x = element_text(size = 17, vjust = -.25)) + labs(color = "")

dev.off()

## Table A.2 ##
m.maj.strong <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
              extremity +
              extremity.lead +
              extremity*extremity.lead +
              (1|state) +
              (1|leg.s),
              data = wr.id[which(wr.id$inparty == 1 & wr.id$maj_agenda_power == 3), ], family = binomial (link = logit))

m.maj.weak <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
              extremity +
              extremity.lead +
              extremity*extremity.lead +
              (1|state) +
              (1|leg.s),
              data = wr.id[which(wr.id$inparty == 1 & wr.id$maj_agenda_power < 3), ], family = binomial (link = logit))

## Table A.3 ##
# Control variables #
# Term limits
wr.id$termlimits <- ifelse(wr.id$state == "me" | wr.id$state == "ca" | wr.id$state == "co" | wr.id$state == "ar" | wr.id$state == "mi" | wr.id$state == "fl" | wr.id$state == "oh" | wr.id$state == "sd" | wr.id$state == "mt" | wr.id$state == "az" | wr.id$state == "mo" | wr.id$state == "ok" | wr.id$state == "ne" | wr.id$state == "la" | wr.id$state == "nv", 1, 0)

# Cosponsorship limits
wr.id$cosp <- ifelse(wr.id$st == "IN" | wr.id$st == "MN" | wr.id$st == "NH" | wr.id$st == "ND" | wr.id$st == "WV" | wr.id$st == "WY", 1, 0)

# Chamber polarization
wr.id$polar <- abs(wr.id$X2010maj_id - wr.id$X2010min_id)

# Four year term dummy
wr.id$fouryear <- ifelse(wr.id$state == "al" | wr.id$state == "la" | wr.id$state == "md" | wr.id$state == "ms" | wr.id$state == "nd", 1, 0)

# Professionalism (http://www.ncsl.org/research/about-state-legislatures/full-and-part-time-legislatures.aspx)
wr.id$fulltime <- ifelse(wr.id$state == "ak" | wr.id$state == "ca" | wr.id$state == "wi" | wr.id$state == "mi" | wr.id$state == "il" | wr.id$state == "oh" | wr.id$state == "fl" | wr.id$state == "pa" | wr.id$state == "ny" | wr.id$state == "ma", 3, NA)
wr.id$fulltime <- ifelse(wr.id$state == "ga" | wr.id$state == "id" | wr.id$state == "ks" | wr.id$state == "me" | wr.id$state == "ms" | wr.id$state == "nv" | wr.id$state=="nm"|wr.id$state=="ri"|wr.id$state=="vt"|wr.id$state=="wv"|wr.id$state=="mt"|wr.id$state=="nh"|wr.id$state=="nd"|wr.id$state=="sd"|wr.id$state=="ut"|wr.id$state=="wy", 1, wr.id$fulltime)
wr.id$fulltime<-ifelse(is.na(wr.id$fulltime)==TRUE, 2, wr.id$fulltime)

# Party size
size <- read.csv("PartySize.csv")
wr.id$MajSize <- numeric(nrow(wr.id))
wr.id$MinSize <- numeric(nrow(wr.id))
for(i in 1:nrow(wr.id)){
wr.id$MajSize[i] <- size$MajSize[which(as.numeric(size$State) == as.numeric(wr.id$st[i]))]
wr.id$MinSize[i] <- size$MinSize[which(as.numeric(size$State) == as.numeric(wr.id$st[i]))]
}
wr.id$PartyAd <- (wr.id$MajSize - wr.id$MinSize)/(wr.id$MajSize + wr.id$MinSize)

# Models #
m.maj.c <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 1), ], family = binomial (link = logit),
             control = glmerControl(optimizer = "bobyqa"))

m.min.c <- glmer(matrix(c(Bills.Waffled, Sponsored.Votes - Bills.Waffled), ncol = 2) ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd +             
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 0), ], family = binomial (link = logit),
             control = glmerControl(optimizer = "bobyqa"))

## Table A.4 ##
# Cosponsored #
m.maj.cosp <- glmer(Bills.Sponsored ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 1), ], family = poisson,
             control = glmerControl(optimizer = "bobyqa"))
             
m.min.cosp <- glmer(Bills.Sponsored ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 0), ], family = poisson,
             control = glmerControl(optimizer = "bobyqa"))

# Bills receiving votes #
m.maj.brv <- glmer(Sponsored.Votes ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd + scale(Bills.Sponsored) +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 1), ], family = poisson,
             control = glmerControl(optimizer = "bobyqa"))
             
m.min.brv <- glmer(Sponsored.Votes ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd + scale(Bills.Sponsored) +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 0), ], family = poisson,
             control = glmerControl(optimizer = "bobyqa"))

# Waffle counts #             
m.maj.wc <- glmer(Bills.Waffled ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd + scale(Sponsored.Votes) +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 1), ], family = poisson,
             control = glmerControl(optimizer = "bobyqa"))
             
m.min.wc <- glmer(Bills.Waffled ~ 
             extremity + 
             extremity.lead +
             extremity*extremity.lead +
             termlimits + cosp + polar + fouryear + fulltime + PartyAd + scale(Sponsored.Votes) +
             (1|state) +
             (1|leg.s),
             data = wr.id[which(wr.id$inparty == 0), ], family = poisson,
             control = glmerControl(optimizer = "bobyqa"))             

